Tadpole

Go back to the Convergence Rates page.

Let us set some global options for all code chunks in this document.

knitr::opts_chunk$set(
  message = FALSE,    # Disable messages printed by R code chunks
  warning = FALSE,    # Disable warnings printed by R code chunks
  echo = TRUE,        # Show R code within code chunks in output
  include = TRUE,     # Include both R code and its results in output
  eval = TRUE,       # Evaluate R code chunks
  cache = FALSE,       # Enable caching of R code chunks for faster rendering
  fig.align = "center",
  #out.width = "100%",
  retina = 2,
  error = TRUE,
  collapse = FALSE
)
rm(list = ls())
set.seed(1982)

1 Import libraries

# Install R-INLA package
# install.packages("INLA",repos = c(getOption("repos"),INLA ="https://inla.r-inla-download.org/R/testing"), dep = TRUE)
# Update R-INLA package
# inla.upgrade(testing = TRUE)
# Install inlabru package
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# Install rSPDE package
# remotes::install_github("davidbolin/rspde", ref = "devel")
# Install MetricGraph package
# remotes::install_github("davidbolin/metricgraph", ref = "devel")

library(INLA)
library(inlabru)
library(rSPDE)
library(MetricGraph)
library(Matrix)

library(dplyr)
library(plotly)
library(scales)
library(patchwork)
library(tidyr)
library(ggplot2)
library(reshape2)
library(sf)

library(here)
library(rmarkdown)
library(knitr)
library(grateful) # Cite all loaded packages

library(latex2exp)
library(plotrix)

rm(list = ls()) # Clear the workspace
set.seed(1982) # Set seed for reproducibility

2 Tadpole graph

Follow this link for more details.

# Eigenfunctions for the tadpole graph
tadpole.eig <- function(k,graph){
  x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2]) 
  x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2]) 
  
  if(k==0){ 
    f.e1 <- rep(1,length(x1)) 
    f.e2 <- rep(1,length(x2)) 
    f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
    f = list(phi=f1/sqrt(3)) 
    
  } else {
    f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2) 
    f.e2 <- sin(pi*k*x2/2)                  
    
    f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
    
    if((k %% 2)==1){ 
      f = list(phi=f1/sqrt(3)) 
    } else { 
      f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
      f.e2 <- cos(pi*k*x2/2)
      f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1]) 
      f <- list(phi=f1,psi=f2/sqrt(3/2))
    }
  }
  
  return(f)
}


# Function to compute the true covariance matrix
gets_true_cov_mat <- function(graph, kappa, tau, alpha, n.overkill){
  Sigma.kl <- matrix(0,nrow = dim(graph$mesh$V)[1],ncol = dim(graph$mesh$V)[1])
  for(i in 0:n.overkill){
    phi <- tadpole.eig(i,graph)$phi
    Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*phi%*%t(phi)
    if(i>0 && (i %% 2)==0){ 
      psi <- tadpole.eig(i,graph)$psi
      Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*psi%*%t(psi)
    }
    
  }
  Sigma.kl <- Sigma.kl/tau^2
  return(Sigma.kl)
}

# Define the graph
edge1 <- rbind(c(0,0), c(1,0))
theta <- seq(from = -pi, to = pi,length.out = 100)
edge2 <- cbind(1+1/pi+cos(theta)/pi, sin(theta)/pi)
edges <- list(edge1, edge2)
graph <- metric_graph$new(edges = edges)

# parameters
h.ok <- 2^-10
type <- "covariance"
type_rational_approximation = "brasil"
rho <- 0.5
#m = 4
sigma <- 1
n.overkill = 1000

# Mesh sizes
h_aux <- seq(5.5, 4, by = -1/2)
h_vector <- 2^-h_aux
h_label <- paste0("2^-", h_aux, "")
h_label_latex <- sprintf("$2^{-%f}$", h_aux)

# Beta values
beta_aux <- c(4, 6, 6.5, 7, 7.5)
beta_vector <- beta_aux/8
beta_label <- paste0(beta_aux, "/8")
theoretical_rate <- pmin(4*beta_vector-1/2,2)

graph.ok <- graph$clone()
# Build graph with overkill mesh
graph.ok$build_mesh(h = h.ok)
graph.ok$compute_fem()

# Get the overkill mesh locations
loc.ok <- graph.ok$mesh$VtE # or graph.ok$get_mesh_locations()

# Initialize the list of graphs and the list of projection matrices
graphs <- list()
A <- list()
for(i in 1:length(h_vector)){
  graphs[[i]] <- graph$clone()
  graphs[[i]]$build_mesh(h = h_vector[i])
  A[[i]] <- graphs[[i]]$fem_basis(loc.ok)
}

cov.error.folded <- matrix(NA, nrow = length(h_vector), ncol = length(beta_vector))

m_values <- c()
for (j in 1:length(beta_vector)) {
  beta <- beta_vector[j]
  alpha <- 2*beta
  fract2beta <- alpha - floor(alpha)
  nu <- alpha - 0.5
  kappa <- sqrt(8*nu)/rho
  tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)))  #sigma = 1, d = 1

  Sigma.folded <- gets_true_cov_mat(graph = graph.ok,
                               kappa = kappa,
                               tau = tau,
                               alpha = alpha,
                               n.overkill = n.overkill)

  for (i in 1:length(h_vector)) {
    h <- h_vector[i]
    m <- min(20, ceiling((min(4*beta - 1/2,2) + 1/2)^2*log(h)^2/(4*pi^2*fract2beta)))
    m_values <- c(m_values, m)

    Sigma <- matern.operators(alpha = alpha,
                             kappa = kappa,
                             tau = tau,
                             m = m,
                             graph = graphs[[i]],
                             type = type,
                             type_rational_approximation = type_rational_approximation)$covariance_mesh()

    Sigma.approx <- A[[i]]%*%Sigma%*%t(A[[i]])

    cov.error.folded[i,j] <- sqrt(as.double(t(graph.ok$mesh$weights)%*%(Sigma.folded - Sigma.approx)^2%*%graph.ok$mesh$weights))
  }
}

print(m_values)

##  [1] 20 20 20 20  5  4  4  3  4  4  3  2  4  3  3  2  3  3  2  2

slope_tadpole <- numeric(length(beta_vector))

for (u in 1:length(beta_vector)) {
  slope_tadpole[u] <- round(coef(lm(log(cov.error.folded[,u]) ~ log(h_vector)))[2], 7)
}

transposed_df <- data.frame(t(data.frame(beta = beta_label, theoretical_rate = theoretical_rate, slope = slope_tadpole)))
rownames(transposed_df) <- c("beta", "Theoretical rates", "Tadpole graph")
colnames(transposed_df) <- NULL
# Display the transposed data frame
transposed_df |> paged_table()

loglog_line_equation <- function(x1, y1, slope) {
  b <- log10(y1 / (x1 ^ slope))

  function(x) {
    (x ^ slope) * (10 ^ b)
  }
}

guiding_lines <- matrix(NA, nrow = length(h_vector), ncol = length(beta_vector))
for (j in 1:length(beta_vector)) {
  guiding_lines_aux <- matrix(NA, nrow = length(h_vector), ncol = length(h_vector))
  for(k in 1:length(h_vector)){
    point_x1 <- h_vector[k]
    point_y1 <- cov.error.folded[k, j]
    slope <- theoretical_rate[j]
    line <- loglog_line_equation(x1 = point_x1, y1 = point_y1, slope = slope)
    guiding_lines_aux[,k] <- line(h_vector)
  }
  guiding_lines[,j] <- apply(guiding_lines_aux, 1, mean)
}

# Generate default ggplot2 colors
default_colors <- scales::hue_pal()(ncol(guiding_lines))

# Create the plot_lines list with different colors for each line
plot_lines <- lapply(1:ncol(guiding_lines), function(i) {
  geom_line(data = data.frame(x = h_vector, y = guiding_lines[, i]),
            aes(x = x, y = y), color = default_colors[i], linetype = "dashed", show.legend = FALSE)
})

df <- as.data.frame(cbind(h_vector, cov.error.folded))
colnames(df) <- c("h_vector", beta_label)
df_melted <- melt(df, id.vars = "h_vector", variable.name = "column", value.name = "value")

p3 <- ggplot() +
  geom_line(data = df_melted, aes(x = h_vector, y = value, color = column)) +
  geom_point(data = df_melted, aes(x = h_vector, y = value, color = column)) +
  plot_lines +
  labs(title = "Tadpole graph",
    x = TeX("$\\hat{h}$"),
       y = "Covariance Error",
       color = expression(beta)) +
  scale_x_log10(breaks = h_vector, labels = TeX(h_label_latex)) +
  scale_y_log10() +
  theme_minimal() +
  theme(text = element_text(family = "Palatino"))

p3